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| Abstract. We study the rate at which stars spiral into a massive black hole (MBH) due to 

the emission of gravitational waves (GWs), as a function of the mass M. of the MBH. In 
the context of our model, it is shown analytically that the rate approximately depends on the 
MBH mass as M^ 1 ^. Numerical simulations confirm this result, and show that for all MBH 
masses, the event rate is highest for stellar black holes, followed by white dwarfs, and lowest 
for neutron stars. The Laser Interferometer Space Antenna (LISA) is expected to see hundreds 
of these extreme mass ratio inspirals per year. Since the event rate derived here formally 
diverges as M, — > 0, the model presented here cannot hold for MBHs of masses that are too 
low, and we discuss what the limitations of the model are. 
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! 1. Introduction 

VD 

The inspiral of a compact remnant into a massive black hole (MBH) leads to the emission of 
gravitational waves (GWs) with frequencies that will be detectable by the Laser Interferometer 
Space Antenna (LISA). Such extreme mass ratio inspirals (EMRIs) will lead to a wealth 
of astrophysical information, such as highly precise measurements of the mass and spin of 
the MBH, the mass of the inspiraling object, the distance to the source, and the direction 
of the source on the sky 0. LISA will be sensitive to MBHs with masses in the range 
10 4 M Q £ M, ^ 1O 7 M , and is expected to detect thousands of events during its mission 
lifetime [14]. The existence of MBHs with masses M, < 10 6 M Q is still hypothetical, in 
spite of considerable circumstantial evidence (see ll37l for a review). Some of this evidence 
stems from ultra-luminous X-ray sources (ULXs). These are sources of X-rays that are more 
luminous than the Eddington luminosity for a 1OM stellar black hole (BH), but are located 
away from the center of their host galaxy, which makes it improbable that they are very 
massive, as dynamical friction would have brought them to the center. Arguments that selected 
ULXs are more massive than stellar BHs include mass estimates from the mass-temperature 
relation ll35l l36l and quasi-periodic oscillations ifTOl l33l . Dynamical analysis also indicates 
that some globular clusters may contain quiescent intermediate mass BHs [fT9j[T6l[39l . 

The rate at which EMRIs form near a given MBH is still rather uncertain, in part due to 
the scope of different dynamical processes that need to be taken into account. Estimates of 
event rates []2S|45l[M[ID[3l[Illl|251|2il27l have typically focused on MBHs with masses 
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of M, ~ 1O 6 M (see 11281 for a review). One reason for this is that the signal for such a MBH 
is strongest, but another, more practical reason, is that in that case the system can be scaled 
with the situation in the Galactic center, which contains a MBH of M, « (3 — 4) x 1O 6 M 
11201 l44l [H |2TJ . The Galactic center has been studied in great detail (see [Q]| for a review), 
and the knowledge we have of the stellar properties of the Galactic center can be exploited to 
calibrate models of stellar galactic nuclei. 

A possible Ansatz to estimate event rates for lower mass MBHs is to assume that the 
empirical M, — a relation 

M. = 1O 8 M (—£ -V (1) 



V 200km s" 1 . 

between the M. and the velocity dispersion of the host galaxy |[9l [T71 |46l extends to lower 
mass MBHs. By making this assumption, inspiral rates for these hypothetical MBHs can 
be estimated. It should be stressed that the continuation of the M, — a relation to lower 
masses has no empirical foundation - indeed, the very existence of lower mass MBHs is not 
universally accepted. 

One obvious interest in the detection of the inspiral into a lower mass MBH by LISA 
is that it would immediately provide unambiguous proof for the existence of such MBHs. 
In addition, an interesting possibility that is unique for lower mass MBH was suggested 

by ll34l|Sesana et al.(2008)Sesana. Vecchio. Eracleous. & Sigurdssonf. Since the tidal radius 

173 

r t = (M./M*) ' .R* of a star with mass M* and radius has a weaker dependence on M, 

than the Schwarzschild radius rs = 2GM,/ c 2 , tidal effects become more important compared 

to general relativistic effects as M, decreases. As a result, when MBHs of mass M. <, 10 5 

capture a white dwarf (WD), the WD may be disrupted by the tidal field during the last phases 

of inspiral. This may lead to a very interesting coincidence of a LISA detection with an 

electromagnetic counterpart. Lower mass MBHs also contribute differently to the stochastic 

background of GWs from EMRIs, due to the GWs emitted prior to the phase where they reach 

a signal to noise ratio high enough for detection [6l . 

These questions motivate us to revisit the dependence of the inspiral rate T on M.. [|25l 

used analytical arguments to find T oc M" 1 / 4 . Here we rephrase and extend these arguments 

for systems that include mass-segregation ($2]). We then perform a series of numerical 

simulations, to derive the dependence T(M.) numerically (©, and show that the dependence 

is in good agreement with the analytical result. We discuss the validity, extensions, and 

consequences of our model in ($4]). 



2. EMRI rate model and scaling with mass 

In order to compensate for our ignorance about lower mass MBH systems, we make the 
Ansatz discussed in the introduction of assuming that the M, — a relation can be extended 
to M. < 10 6 M Q , and we assume that the mass enclosed within the radius of influence 
r h = GM./a 2 is proportional (and comparable) to that of the MBH. Based on these 
assumptions, we derive an analytical expression (fl~2l) for the dependence of the event rate 
on M, and present a numerical model to calculate these event rates. 
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2.1. Scaling galactic nuclei with mass 

With the aforementioned assumptions, we find a number of scaling relations. The radius of 
influence, where the MBH dominates the potential, is given by 



GM, 



a 

( M V' 2 

The fact that the enclosed mass of stars within r h is proportional to the MBH mass 
implies that the number density at the radius depends on M, as 



n, = 4xl0 4 pc" 3 ( 

(for the normalization, see lfi~8l0 . 

The relaxation time then increases with the mass of the MBH as 



3(2trx 2 ) 3 / 2 
h Z2n 2 G 2 M 2 n h \xiK 

( M \ 5/4 

= 6G *<(t7Wm;) ■ (4) 

In this expression, the Coulomb logarithm is approximately In A = In (M./M*); for our 
analytical arguments we discard the weak mass -dependence of the Coulomb logarithm. 



2.2. Model for EMRI rate 

We follow the analysis presented in E71 . who solve the Fokker- Planck equations for a multi- 
mass model with 4 species: main sequence stars (MSs) of mass Mms = M Q ; WDs of mass 
Mwd = O.6M ; neutron stars (NSs) of mass Mns = 1.4M ; and stellar BHs of mass 
^bh = IOMq. In all cases we assume that the stars far away from the MBH have equal 
velocity dispersions a, and that the unbound stars have a distribution function 

g M {x) = C M exp(x) (x < 0), (5) 

where x = E/a 2 is the dimensionless energy of the stars, and E is the negative specific energy 
with respect to the MBH. The prefactor Cm determines the fraction of compact remnants of 
type M in the unbound stellar population. The distribution © is what may be expected if the 
system was formed by violent relaxation; the stellar distribution of stars bound to the MBH 
is not very strongly affected by the exact form of this distribution (although it is dependent 
on its normalization), see e.g. 0. To simplify our analysis, we do not consider the evolution 
of unbound stars, and assume that they are in steady state. This is mainly a simplification, 
but there is some justification in the fact that the relaxation time far away from the MBH 
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increases rapidly with distance, so that the dynamical evolution of these stars is slow. Also, 
there is likely star formation far away from the MBH, which will be harder to model, so that 
it is not clear how much realism would be obtained by including dynamics far away from 
the MBH. We note that the stellar models presented in ifTTl [T3l do include dynamics for all 
stars in the system. For our models, we will assume that Cm = (1, 0.1, 0.01, 0.001) for MSs, 
WDs, NSs and BHs, consistent with models of continuous star formation that appear to be 
consistent with Galactic center observations dHH]. 

The Fokker-Planck equations include a loss-cone term, and in dimensionless units, the 
equations have the following form (for details, see E71l : 1131 ) 

dg M (x,T) 5/2 d 

dr = ~ X dx 1 ~ ^ 

where Rm{x) is the rate at which stars of type M, with energies in the interval (x, x + dx), 
reach the loss-cone and are eaten by the MBH. The rate at which stars flow through a given 
energy x is given by Q = Q[g(x, r),x); this non-linear term represents 2-body relaxation. We 
stress that these equations do not depend on the MBH mass (apart from the assumption that 
it is much more massive than any of the individual stars). The same Fokker-Planck equation 
therefore applies to all systems, regardless of their MBH. The mass does enter, of course, 
in the dimension-full equations, and dictates for example how many years it takes to reach 
steady state. 

The rate at which stars enter the loss-cone is much higher than the EMRI formation 
rate: many stars plunge prematurely without spiraling in. These stars never reach an orbit 
detectable by LISAfl 

In dimension-full units, the rate at which stars with semi-major axes less than a come 
into the loss-cone is then 

T M (< a) = 2V2I Q ( a/r ' h dyy l ' 2 R M (y) , (7) 



o 



where the flow I is given by 




_ 8tt 2 3 (GM±) 2 \nAn h 

J = ^lK r h n h 1 

a 6 

-1/4 

(8) 

which is approximately the number of the stars within r h divided by the relaxation time. The 
actual flow is in some situations much smaller than I (i.e., Jq 1 dyy l / 2 R,M (y) «C 1), as 
was realized by 0, who derived several analytical results from such "zero-flow" solutions. 
However, it was shown by that in cases such as here, where there are very rare, massive 
objects (i.e., stellar BHs), maximal flows can in fact be achieved. 

For LISA, the only objects that are accreted by the MBH and detected are those 
(compact) stars that spiral into the MBH gradually, without plunging prematurely into the 

| An exception to this is the Galactic center where gravitational bursts from fiy-bys are possibly detectable 

i@a. 
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loss-cone. Therefore, for in-spirals, the integrand in equation C7]) should be convolved with 
a function Sm{cl) that gives the probability for a star to spiral in without plunging. This 
function filters out all stars that plunge into the MBH without leading to a detectable EMRI 
source ("plunges"), and selects only the stars that do lead to a slowly inspiraling star. We 
compute this function with the Monte Carlo methods as described in 11251 (see also ll23l . who 
used a very similar method). The total EMRI rate is then given by 

1 m{< a) = 2V2I 1 dyy 1 ^S M (y)R M (y). (9) 
Jo 

The function Sjvf (a/r^) selects only those captured stars that spiral in successfully, which 
is a small fraction of the total number of captures [so J^ Th dyy 1 ^ 2 SM(y)RM (y) "C 
Jt h dyy 1/2 R M (y)l of order 10% (2). 

In contrast to the solution of the dimensionless Fokker-Planck equations, the function 
Sm(cl) may depend on the MBH mass, because the loss-cone, which is determined by the 
angular momentum of the last stable orbit Jlso — 4GM./c, depends on M„ and provide an 
additional length scale to the problem. However, it was argued by [|25l , that Sm{cl) can be 
approximated by a step-function that is zero for a > a c , with 

3/(3-2p M ) 



n 



(10) 



where 



, / SVGM^TA 2 ^ 85tt M,c 2 

d c = I - 2 ; E 1 = 77—^ aj , (11) 



TCC 



3x2 13 M. 



and pm is the slope of the distribution function ^m(^)- It follows that d c oc M^ 2 , and hence 
that a c /rh is independent of M, (see also figured)). Hence, at the level of this approximation, 
it follows that the integral in equation © does not depend on mass, and thus that 



T M = 2V2I / dyy 1/2 S M (y)RM (y) 
Jo 

w 2V2I / dyy^Ru (y) 
Jo 

oc J oc M" 1/4 . (12) 

The result (fT2l) is derived from the assumption that S]\i(y) may be approximated by a 
step function Sm(v) ~ 9 ( a c/rh — V) with a c /r^ as in equation (flOl) . This result is based 
on several approximations compared to the numerical Fokker-Planck / Monte Carlo model. 
To test the validity of this approximate expression, we perform Monte Carlo experiments to 
determine the function Sm- 

The calculation for the EMRI rate for all systems we consider here then goes through the 
following steps, (i) calculate the steady state of equations © once for all systems, (ii) From 
the outcome, find the dimensionless relaxation time 
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(iii) Convert this dimensionless time-scale into a dimension-full time through equation ©, 
and use this time in the Monte Carlo simulations to find the function Sm{o)- (iv) Calculate 
the dimension-full integral © to calculate the EMRI rate. 

3. Numerical results 

We consider MBH masses of M. = 3 x 1O 4 M , M, = 3 x 1O 5 M , and M, = 3 x 1O 6 M , 
with stellar populations as indicated in $2l 

In figure (Q~|), we show the inspiral fractions Swd(o) for WDs for the MBH masses 
considered. Inspiral probabilities are high for small semi-major axes, as anticipated. 
The analytical arguments given in the previous section suggest that Swd^AVi) does not 
strongly depend on the MBH mass. To test this idea, we multiply the horizontal axis by 
y^3 x 1O 6 M /M. in order to extract the MBH dependence of r h (equation [2]). The curves 
now coincide approximately, confirming that the dependence on M. is weak. The figure also 
shows that 5 , WD (a) can be reasonably approximated by a step-function. 
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Figure 1. Inspiral fraction 5wd(«) for WDs, for several MBH masses. The solid lines show 
the direct results, whereas the crossed data points were re-scaled by y3x 10 6 M Q /M.. The 
inspiral rate is high for small semi-major axes, and drops to zero at large semi-major axes, with 
<Swd = 0.5 at ~ 2mpc (M./3 x 10 6 M Q ) . The fact that the scaled lines coincide with 
the line for M, = 3 x 10 6 Mq implies that S]\i(a/rf l ) is to good approximation independent 
on the MBH mass. 

The resulting event rates are reported for all different species in table (OQ). The numbers 
for M, = 3 x 10 6 are different by a factor £ 2 compared to those in 11271 . probably mainly 
due to slightly different boundary conditions. 
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Table 1. EMRI rates 



M 


WD 


NS 


BH 


[M e ] 


[Gyr- 1 ] 


[Gyr- 1 ] 


[Gyr- 1 ] 


3 x 10 4 


62 


26 


710 


3 x 10 5 


34 


13 


650 


3 x 10 6 


20 


7 


400 



The event rates are plotted for the different species as a function of mass in figure ©. 
This figure shows that the event rate indeed decreases with M„ and the relation T M oc M, 1 ^ 
is confirmed. 




1 1 ■ — ■ — ■ — 1 

10000 100000 le+06 le+07 

M[M ] 

Figure 2. EMRI event rates as a function of MBH mass for BHs (top line), WDs (middle line), 
and NSs (bottom line). The straight lines give approximate power laws. For WDs and NSs the 
slope of —0.25 is in very good agreement with the analytical result, while for BHs, the slope 
of —0.15 is slightly flatter. 



4. Summary and discussion 

We analyzed the EMRI rate for several MBH masses in the context of the 11251 and E71 
models. The rate depends on MBH mass as oc M~ 1//4 . The analytical mass dependence 
(equation [121 is confirmed by Monte Carlo simulations. 

It is not trivial to calculate the projected EMRI detection rate from a given event rate 
lfl4"l . The inspiral rate estimated here would give rise to a detection rate of ~ 10 3 over the 
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LISA life time of 5 years IfBH . 

Since the event rate increases as the M. mass decreases, the assumptions for the models 
here must break down at some point, since otherwise the MBH would accrete more than its 
own mass within the age of the system. This is already the case for the M, = 3 x 1O 4 M 
considered here. Several possible solutions to this paradox are that (i) such MBHs do not 
exist, or exist only during transient phases much shorter than the age of the Universe (for 
example, because of such high capture rates); or that (ii) somehow the assumed boundary 
conditions are not fulfilled, for example because the M, — a relation is not satisfied MBHs 
with M, < 1O 6 M , or BHs do not sink in effectively from large distances, and are not formed 
by new star formation, so that the boundary condition for the unbound stars is not maintained 
over time. In the latter situation, the WD and NS rates, which are lower (and bring in less 
mass), could still be correct, and perhaps the actual rates would be even higher if stellar 
BHs do not push these stars out by mass segregation. To test this idea, we performed a 
simulation with a M. = 3 x 1O 4 M MBH, similar to the ones presented in the previous 
sections, but without any stellar BHs. In this example, the event rate increases by a factor 
~ 2.5 to T WD = 146Gyr _1 . We finally note that we find that the MBH grows as M, oc t 4/5 
as a result of EMRI captures, much slower than the exponential growth that would result from 
Eddington limited accretion. 

In a recent paper, SOI study the stellar BH distribution around MBHs in a different 
context, where they consider cases with a much flatter mass function than we assume here, 
and with more massive stellar BHs. This will likely increase the EMRI rate further (and hence 
restrict further the MBH range for which the analysis is valid). Here we do not consider this 
situation. 

The models presented here do not include the effects of resonant relaxation ||42]|. 
Resonant relaxation increases the rate at which the angular momenta of the stars randomize, 
but the effect is largely destroyed for very eccentric orbits due to general relativistic 
precession. It was argued by [26] that resonant relaxation will likely increase the EMRI 
rate by a factor of £10. Many aspects of resonant relaxation remain unclear, however. 
Although the eccentricity dependence of random torques in a power-law cusp was found by 
ll22l . the coherence time of the torque, and its dependence on eccentricity, remain unknown. 
In addition, there is the related effect of resonant friction fl42|. which has not been considered 
so far in the literature in detail. It can be concluded that resonant relaxation will probable 
modify the EMRI event rate, but a more detailed analysis is premature at this point. 

It has been suggested that triaxiality can also affect the EMRI rate |[24|. Angular 
momentum is not a conserved quantity in a triaxial potential, and as a result it can potentially 
change much faster than can be achieved by relaxation mechanisms. However, it seems 
unlikely that the potential is significantly triaxial within the radius of influence, since the 
potential is there dominated by the MBH. Since Sm{o) = for a > 10~ 2 pc (see figure [B, 
enhanced evolution of angular momenta outside this distance would not lead to any additional 
EMRIs: if more stars are captured by the MBH, these captures will all result in plunges. 
Triaxiality outside the radius of influence can therefore only increase the capture rate of 
plunges. Similar conclusions apply to massive perturbers [|4Tll48l . 
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We have focused on the direct capture of compact remnants by MBHs. An alternative 
route to making a LISA source with an intermediate mass black hole is by tidal capture of a 
main sequence star which circularizes near the MBH, and after it leaves the main sequence 
spirals into the MBH as a compact remnant [[3T1I301 . Such different formation mechanisms 
may lead to higher event rates. 
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